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We compare eigenvalue densities of Wigner random matrices whose elements are independent 
identically distributed (iid) random numbers with a Levy distribution and maximally random ma- 
trices with a rotationally invariant measure exhibiting a power law spectrum given by stable laws 
of free random variables. We compute the eigenvalue density of Wigner-Levy (WL) matrices using 
(and correcting) the method by Bouchaud and Cizeau (BC), and of free random Levy (FRL) ro- 
tationally invariant matrices by adapting results of free probability calculus. We compare the two 
types of eigenvalue spectra. Both ensembles are spectrally stable with respect to the matrix addi- 
tion. The discussed ensemble of FRL matrices is maximally random in the sense that it maximizes 
Shannon's entropy. We find a perfect agreement between the numerically sampled spectra and the 
analytical results already for matrices of dimension TV = 100. The numerical spectra show very 
weak dependence on the matrix size TV as can be noticed by comparing spectra for TV = 400. After a 
pertinent rescaling spectra of Wigner-Levy matrices and of symmetric FRL matrices have the same 
tail behavior. As we discuss towards the end of the paper the correlations of large eigenvalues in the 
two ensembles are however different. We illustrate the relation between the two types of stability 
and show that the addition of many randomly rotated Wigner-Levy matrices leads by a matrix 
central limit theorem to FRL spectra, providing an explicit realization of the maximal randomness 
principle. 

PACS numbers: 02.50.-r, 02.60.-x, 89.90.+n 



I. INTRODUCTION 

Applications of random matrix theory cover many branches of physics and cross-disciplinary fields [l[ involving 
multivariate analysis of large and noisy data sets |2|. The standard random matrix formulation belongs to the 
Gaussian basin, with a measure that is Gaussian or polynomial with finite second moment. The ensuing macroscopic 
spectral distribution is localized with finite supports on the real axis. The canonical distribution for a Gaussian 
measure is Wigner's semi-circle. 

The class of stable (Levy) distributions [|| is however much larger (the Gaussian class represents only one fixed 
point in the stability basin of the Levy class), and one is tempted to ask, why the theory of random Levy matrices is 
not so well established. The case of Levy randomness is far from being academic, and many distributions in physics 
and outside (finance, networks) exhibit power- like behavior referred to as fat or heavy tails [J]. 

One of the chief reasons for why the theory of random Levy matrices is not yet well understood is the technical 
difficulty inherent to these distributions. First, even for one-dimensional stable distributions, the explicit form of the 
probability distribution functions (pdf) is known analytically only in few cases [1, @]. Second, Levy distributions have 
divergent moments, and the finiteness of the second moment (condition for the Gaussian stability class) is usually a 
key for many of the techniques established in Random Matrix Theory. Third, numerical studies involving power-like 
behavior on infinite supports require enormous statistics and are very sensitive to systematic errors. 

In 1994, Bouchaud and Cizeau 0] (hereafter BC) considered large TV x TV random symmetric matrices with entries 
sampled from one-dimensional, stable distributions. In the large TV limit they obtained analytical equations for the 
entries of the resolvent, and then checked their predictions for the spectra using numerically generated spectra by 
random sampling. The agreement was fair, although not perfect. Contrary to the standard Gaussian-like ensembles, 
the measure in the BC approach was not rotationally invariant. 

In 2002, following the work in Q, we have suggested another Levy- type ensemble @ (hereafter FRL). By construc- 
tion its measure is rotationally invariant. The average spectral distribution in this ensemble is stable under matrix the 
convolution of two independent but identical ensembles. It is similar to the stability property of one-dimensional Levy 
distributions. The measure is non-analytic in the matrix H and universal at large H with a potential V(H) sa ln_ff 2 . 
This weak logarithmic rise in the asymptotic potential is at the origin of the long tail in the eigenvalue spectra. 

The present work will compare the WL and FRL results as advertised in 10]. In section 2 we reanalyze and correct 
the original arguments for the resolvent presented in 0] • Our integral equations for the resolvent and spectral density 
are different from the ones in Q ■ We carry explicit analytical transformations and expansions to provide insights to 
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the spectrum. We show that there is a perfect agreement between the analytical and numerical results obtained by 
sampling large Levy matrices. We also discuss the relation of ours and BC's results. In section 3 we recall the key 
concepts behind FRL ensembles. In large N the resolvent obeys a simple analytic equation. The resulting spectra are 
compared to the spectra following from the corrected BC analysis. The WL and FRL matrices represent two types 
of stability under matrix convolution. In both cases we have a power behavior in the tails of the spectrum. By a 
pertinent rescaling we may in fact enforce the same tail behavior and compare the spectra. The observed differences 
disappear in the Gauss limit and become more pronounced in the Cauchy limit. In sectin 4 we explain the relation 
between the two types of stability on a simple example: we construct sums of WL matrices rotated by random O(N) 
matrices and show that the spectrum of these sums converges by a matrix central limit theorem to the pertinent 
symmetric FRL spectrum. Our conclusions are in section 5. 



II. WIGNER-LEVY MATRICES 
Definition of the ensemble 

In a pioneering study on random Levy matrices, Bouchaud and Cizeau Q discussed a Wigner ensemble of N x N 
real symmetric random matrices with elements being iid random variables: with probability density function following 
a Levy distribution P{x) = N 1 ^^^ {N x ^x), with /x being the stability index, (3 - the asymmetry parameter, and 
C - the range of the distribution (see below). We shall call these matrices Wigner-Levy (WL) or Bouchaud-Cizeau 
(BC) matrices. The probability measure for the ensemble of such matrices is given by: 

dn WL {H) = Y\P{H l3 )dH l3 (1) 

i<j 

The scaling factor N 1 ^ in pdf makes the limiting eigenvalue density independent of the matrix size N when N — > oo. 
Alternatively one can think of the matrix elements as if they were calculated as = hij/N 1 ^ with hij being 
iid random numbers independent of N: p(x) = L^'^(x). 

Levy distributions are notoriously hard to write explicitly (except in few cases), but their characteristic functions 
are more user friendly [f| 

L c ^{x) = ±- J dkL(ky k * (2) 

where the characteristic function is given by 

logL(fc) = -C\k\^(l + i/3 sign(fc) tan(7r/i/2)). (3) 
The parameters fi, (3 and C are related to the asymptotic behavior of L^{x) 

(4) 

with the /i-dependent parameter 7(/i) given by 

7 ( M )=r(l+ M )sin(^) (5) 

Here /i is the stability index defined in the interval (0, 2], — 1 < f3 < 1 measures the asymmetry of the distribution 
and the range C > is the analogue of the variance, in a sense that a typical value of x is C 1 ^. A standard choice 
corresponds to C = 1. 

We shall consider here only the stability index in the range (1,2), although as will be shown later, results obtained 
in this range seem to be valid also for fi = 1 . We also assume that all random variables have zero mean. 



Determination of the eigenvalue density 



A method of calculating the eigenvalue density of the Wigner-Levy matrices was invented by Bouchaud and Cizeau 
7]. Let us in this section briefly recall the main steps of the method. It is convenient to introduce the resolvent, 
called also Green's function: 

g(z) = ±{TrG(z)) (6) 
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where elements of the matrix G(z) are 

G ij (z) = (z-H)^ (7) 

and the averaging is carried out using the measure |T]). The resolvent contains the same information as the eigenvalue 
density p(X). Indeed if one approaches the real axis one finds that p(A) = — l/7rlim e ^ + Im g(X + ie). This is how 
one usually calculates p(X) from g(z). For the Wigner-Levy ensemble, where individual matrix elements have large 
scale-free statistical fluctuations a slightly different method turns out to be more practical - a method which allows 
one to avoid problems in taking the double limit (first N — > oo and e — ► + ) in which the fluctuations are suppressed 
in an uncontrollable way in the presence of an imaginary part of z. In the BC method z is kept on the real axis and 
fluctuations are not suppressed, so one can safely take the large N limit. 

The method goes as follows One first generates a symmetric N x N random matrix H using the measure ([I]) 
and then by inverting H— z one calculates the resolvent G(z) j7]). Next one adds a new row (and a symmetric column) 
of independent numbers identically distributed as those in the old matrix H. One obtains a new (7V + 1) x (N + l) 
matrix H + , where +1 emphasizes that it has one more row and one more column than H. It is convenient to number 
elements of the original matrix Hij by indices running over the range i,j — 1, . . . , N, and assign the index to the 
new row and the new column so that now indices of H +l run over the range 0, . . . , N . If one inverts the matrix 
H +1 — z one obtains a new (N+l) x {N+l) resolvent G +1 (z). One can show that it obeys a recursive relation [jj 

Z ~chz) =Hm+ ^ H ^ G ^ = wh + t ^§^- (8) 

which relates the element G c J 1 (z) of the [N+l) x (N+l) resolvent to the elements of the old N x N resolvent G(z). 
One can also derive similar equations for off-diagonal elements of G +1 (z): 

G tj( z ) _ hojGjj(z) 
Gti(z) V NV " 

The difference between the probability distribution of the elements of G +1 (z) and of G(z) dissapears in the limit 
N — > oo. The diagonal elements of the matrix G +1 (z) are identically distributed as the diagonal elements of the 
matrix G(z). The same holds for off-diagonal ones. Moreover in this limit all elements of the G matrix become 
independent of each other. In particular all diagonal elements of G(z) become independent identically distributed 
(iid) random variables as N — ► oo. One can use equations © and ^ to derive self-consistency equations for the 
probability density function (pdf ) for diagonal and the pdf for off-diagonal elements. We are here primarily interested 
in the distribution of the diagonal elements, as we shall see below. The self-consistency equation for the probability 
distribution of diagonal elements follows from the equation ([8|) and is independent of the distribution of the off-diagonal 
elements. This can be seen as follows. Let us first define after Bouchaud and Cizeau Q a quantity: 

S a (z) = z-— 1— (10) 

It is merely a convenient change of variables suited to the left hand side of equation ©. It is clear that if one 
determines the pdf for S = Sq(z), one will also be able to determine the pdf for G = Gqq(z) since the two pdfs can 
be obtained from each other by the change of variables (fT0|) : 



where Pg and Ps are pdfs for Gqq(z) and Sq(z) respectively. Since the probability distribution Pq is identical for all 
diagonal elements of G, it remains to determine the probability distribution Ps for the quantity Sq(z). 

Let us sketch how to do that. First observe that the first term in the equation ([5]) can be neglected at large N, so 
the equation assumes the form: 

c / \ \ ^ hlfij^z) \ - hoihojGij(z) 

Now observe that by construction hoi are independent of Gij(z). We shall now show that for large N the first term 
in (I12p dominates over the second one. We shall modify here the argument used in [7J, where it was assumed that 
the off-diagonal elements GijN(z), i ^ j are suppressed by a factor 1/N 1 ^ with respect to the diagonal elements. 
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Instead we note that by construction the quantities Gij(z) are statistically independent of hoi, i = 0,...,N. Since 
we shall only be interested by the diagonal elements of the resolvent matrix, we may replace the contribution of the 
off-diagonal elements hoihojGij(z), i ^ j by their averaged values. As a result, the contribution of the off-diagonal 
terms averages out. Note that this is also true in the Gaussian limit \i = 2 where following [?J we may replace all 
elements of (fT2|) by their respective averages. Taking this into account and omitting the subleading contribution Hqq, 
we get 

i 

Thus the problem was simplified to an equation where the left hand side (So = z — 1/Gqq) and the right hand side 
depend only of the diagonal elements of the G-matrix, which as we mentioned before, are identically distributed in 
the limit N — > oo. Using (|13[) one can derive a self-consistency equation for the probability density function (pdf) Pq 
for diagonal elements of the matrix G. 



Generalized central limit theorem 



To proceed further we apply with Bouchaud and Cizeau Q the generalized central limit theorem to derive the 
universal behavior of the sum on the right hand side of (fT3| in the limit N — > oo: 

• i. If the hoi's are sampled from the Levy distribution , the squares ij = for large t, are distributed 
solely along the positive real axis, with a heavy tail distribution: 



(14) 



t { 

irrespective of (3. The sum 

N 



C l 

is distributed following L L (iA. The range parameters C and C are related by (13]) 

2C' 7 W2) = C 7 (/i) (16) 

The factor 2 on the left hand side corresponds to the sum (1 + j3) + (1 — (3) appearing as a contribution from 
positive and negative values of the original distribution of hoi. This relation is important when comparing to the 
numerical results below where C = 1 is used. From now on (and to simplify the equations) we assume instead 
that C = 1. 

• ii. By virtue of the central limit theorem the following sum 

E Gii(z)tj 

i 

of iid heavy tailed numbers ti is for Gu(z) = O(N ) and N — * oo Levy distributed with the pdf: lF^' 13 ^, 
which has the stability index /j/2 and the effective range C(z) and the asymmetry parameter (3(z) calculated 
from the equations: 



N 

X 

and 



1 N 

C(z) = -Y.\G^)r /2 (18) 



, x wJ27 \Gii{zW /2 sign{G u (z)) , 
B(z) = N ^ l ' \ r ' 1 V V " . 19 

*ETlG«(*)K» 



which follow from the composition rules for the tail amplitudes of iid heavy tailed numbers ti defined above. 
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Integral Equations 



We saw in the previous section that the generalized central limit theorem implies for large N that the "self-energy" 
S = Sq(z) is distributed according to the Levy law Ps(S) — L C Jj^'^^ z \s) with the stability index fi/2 being a half 
of the stability index of the Levy law governing the distribution of individual elements of the matrix H, and with 
the effective range parameter C(z) and the asymmetry parameter (3{z) which can be calculated from the equations 
(dHJ) and (|19p. respectively. One should note that the effective parameters C(z),(3(z) of the distribution Ps(S) are 
calculated for C = 1 and that they are independent of (3 of the probability distribution: L^(hij) of the if- matrix 
elements. 

The sums on the right hand side of the two equations (fl8|) . ([T9| for C(z) and j3(z) have a common form 
~k Si f(Gu{z)). Since in the limit N — > oo, the diagonal elements become iid, the sums on can be substituted 
by integrals over the probability density for Gu\ 

1 f;</(G«(z))) = JdG P G (G) }{G) = J ^P S (z I) /(G) (20) 

where in the second step we used the equation (fTTj) . Since the distribution Ps(S) — L C Jj^' ,l3tyZ \S) is known up to 
the values of two effective parameters C(z) and /3(z) the equations (fP8| and (fl9|) can be written as self-consistency 
relations for (3(z) and C(z): 



= f %\G\^L^\z-l/G) 

J —OO 

roc dG\ G \»/2 si {G)L C(z)M*) {z _ 1/G) 
P(Z) = rr " — ■ ( 21 ) 

The symbol f- stands for principal value of the integral. Notice here the difference between our second equation and 
that in Q. In addition to 0] we also note that the resolvent takes the form: 

9(z) = ^J2G u (z) — jf^fjG^f'^W (22) 
The integrals (f2T|) and ([22]) can be rewritten using the new integration variable x — 1/G as 



+00 



C(z) = I dx\x\-»/ 2 L°}?M*\z-x) 



= L + r^ ^(x)\x\-^L°j?>«'\z- X ) 
J+-dx\x\-^L C J^\z-x) 



and 



All the steps above require both z and Gu(z) to be strictly real. The argument cannot be extended to the complex z 
plane. So all integrals above should be interpreted as principal value integrals, wherever it is necessary. Since fi < 2 
all integrals in (j2"5j) are convergent also in the usual sense. 
The equation for g(z) can be rewritten as 



dx 



g(z) = f -=-L C J^ z \x). (25) 



x 



V/2 



Notice the nontrivial dependence on z in the parameters of the Levy distribution. Above equation can be a source of 
confusion, since its structure resembles another representation of g{z) 

poo J\ 

9{z) = f rP(A) (26) 
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which superficially looks as if one could identify in ([25]) x and A and L C ^^ z ^(x) with p(X). This is not the case, and 
one should instead invert (|25l) using the inverse Hilbert transform: 



1 f°° dz 

PW = -{ M. (27) 

In other words, one first has to reconstruct numerically the real part of the resolvent, and only then compute numer- 
ically the spectral function p(X), using the "dispersive relation" (j2"T)) . This is a difficult and rather subtle procedure. 
In the next section we give some analytical insights to the integral equations that would help solve them and extract 
the spectral function. 

Analytical properties useful for numerics 

One cannot do the integrals from previous section analytically. As mentioned, one cannot even write down an 
explicit form of the Levy distribution. It is a great numerical challenge to solve the problem numerically even if 
all the expressions are given. One realizes that already when one tries to compute the Fourier integral ((2)) of the 
characteristic function since one immediately sees that the integrand in the form @ is a strongly oscillating function 
making the numerics unstable. Fortunately using the power of the complex analysis one can change this integral to 
a form which is numerically stable. So in this section we present some analytic tricks which allow one to reduce the 
problem of computing the eigenvalue density as formulated in the previous section to a form which is well suited to 
the numerical computation. 

To simplify (|23[) we proceed in steps. First, we make use of the Levy distribution through its characteristic 

i /*°° 

L CfM = J_ dfce ^ e -C| fc r/ 2 (l + ,Csign(/c))_ (28) 



with £ = (3 tan(7r/x/4). By rescaling through 

k = C~ 2/fl k', (29) 
x = C 2 ^x', 
z = C 2 ^z', 

we can factor out the range L°f 2 {x) = C^^L^x'). Second, we make use of the following integrals, 
dx 



e lkx = -2ie lfcz signfc, 



/ 

" dx , COsjfcV) = |fcr /2-l r(1 _ M/2) Sin(7r/V4); (30) 



o x'^l 2 



[ dx '^fi = |fcT /2 - 1 Bign(fc')r(l-M/2)cos( 7 r M /4). 



Last, we make use of the change of variables p — k ,f1 / 2 . With this in mind, we obtain 

C 2 (z') = i-r (l - f) sin (^) J°° dpcos{p 2 ^z' - COO P)e~ P , (31) 

c( ,j _ JiT dpsm{p 2 ^z' -C(z') p)e~P 
J dpcos(p 2 /^z' — C,(z') p)e~P' 

with C(z') — tan(7r/z/4) (3{z'). For every z 1 we can iteratively solve the equation for ((z r ), then we determine C(z') 
and use z = C 2 ^(z')z' to express everything in terms of z. These transformations solve (|2"3")l . Using the same method 
we rewrite the equation for g{z) as 

9 r°° 

g(z') = C(z') 2 ^g{z) = - dp p( 2 -^/^ sin(p 2 ^z' - p ((z')) e' p . (33) 



M Jo 
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These integral forms are useful to study the small- z' limit. In this case C(z') is an antisymmetric function of z' 
and has an expansion in powers of z' . Using ((z r ) = k^z' + 0(z' 3 ) we can recursively obtain the coefficients of this 
expansion. The first term is k\ = T(l + 2//i)/2. Similarly C(z') is a symmetric function in z' 



c2 ( 2 ') = ^ r ( 1 -f) si °(T) + °^» < 34 > 

For \z'\ large (z' — > ±00) a different approach is needed. In this case we follow Nolan [ll[ and treat the two integrals 
(numerator and denominator) of (l32|) together, i.e. we consider the integral 

dpe~ h ^\ (35) 





where 

h(p) = (1 - i()p + iz'p 2 ^. (36) 

Nolan's idea is to close the contour of integration in the complex p plane in the following way: at p — > 00 we add an 
arc and afterwards continue until p = along the line where Qh(p)=0. Using the parameterization p = re lS we get 
the parametric equation for r(9) along this line, valid for z' > 

(\ ///(2 — /i) 
sin(flo-fl) \ 
z'cos(*o)cos(**)J ■ (3?) 

The angle 9 for the curve we need (p, > 1) is bounded between 9 = arctan(£(z')), where r = and —9\ = —tt/i/A, 
where cos(^-#) is zero. Let us now introduce a new variable ip through 9 = 9q — ip an d < ^ < 6*o +6*1- In this range 
we have 

ycos(6» o )cos(^(-0 - O )) J 



/ 1 \ m/^-m; 



"^ a -")„, „ cos(2^t/. - |6» ) 



008(^0)003(^(^-6)0)) 



After some manipulations we obtain 

l*/(2-n) 



sin(p 2 ^z' - C(z') p)e" p = y # f -J y(^)e _,Wl W (39) 

/ 2 COS (2_A£(^_0 O )) ^ gin0 \ 



r f6 +S 1 / -, \ n/(2-n) 

J cos(p 2 ^z' - C(z') p)e- p = J # I- V(iP)e~ mw 

( 2 sin(2=£(^-0 o )) /l* cos 6^ 



2 — /i cos(-(^> — O ) 2 — /1 sin V' 



I 2-// cos(|(V> - #o)) 2-fisxaipj' 

The resulting integrals look complicated, however they contain both the small-z' and the large-z' asymptotics. 
For z' — > we have "0 = z'p 2 ^ -1 , which reproduces the small-z' expansion presented above. For z' — > 00 we have 
i/> = 6 + 9 X - u/z'^/ 2 . Note that in this limit 

2 2 11 

cos(-O-0 o ))=sm( - " ) (40) 

and the 2' dependence in Re/i(V>) vanishes in leading order. 

The large-z' asymptotics requires some work. For the leading orders we have 
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FIG. 1: Theoretical (black) and numerical (red) eigenvalue distributions for fj, — 1.95 



sin(p 2 ^z' - C(z') p)e- p « r(l + (i/2) sin 9 X {z')~^ 2 
cos(p 2/ ^V - CO') p)e~ p » r(l + ^/2)cos6» 1 (z')-"/ a . 
Thus for z' — > oo we have 



(41) 



and 



C(z') = tan(0i) + 0(1/^ 2 ) 

C(z') = z'-^ /4 (l + C(l/z^ /2 )) 
* = + 0(l/z'"/3)). 



(42) 



(43) 



All the formulas above apply to the case z' > 0. One can also derive similar formulas for z' < 0, it is however more 
practical to use the symmetry properties of the functions C(z') and C( z ')- 

As a final check let us compute g(z'). A rerun of the above transformations on the integrals give 



m = 



( 2 1 



2/(2-n) 



(44) 



/zcos(~(^-0 o )) 



sin(^-^o) 
sin-0 



Notice that both asymptotics follow from this representation. For z' ~ * oo we have g(z') = l/z' + • • • . which implies 
g(z) = l/z + • • • . This can be viewed as a check of the correct normalization of the eigenvalue density distribution. 



Numerical Comparison 

In this section we show perfect agreement between the theoretical analysis for the eigenvalue distribution (|27|) 
and the numerically generated eigenvalue distribution. The latter is obtained by diagonalizing N x N random Levy 
matrices sampled using the measure ([1}. The former was generated by calculating numerically g(z) as detailed above. 
We performed numerically the inverse Cauchy transform (|27[) . It is important to note that the integral transforms 




A 

FIG. 4: Theoretical (black) and numerical (red) eigenvalue distribution for /i — 1.25 
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FIG. 5: Theoretical (black) and numerical (red) eigenvalue distribution for /j, — 1.00 



entering into the definition of the eigenvalue density p(A) in (|27|l converge slowly for z — > oo. For that, we have used 
the asymptotic expansion of g(z) to perform the large- z part of the integrals. 

As noted above, all the above analytical results were obtained using a specific choice of the scale factor (|T5|) . The 
comparison with the numerically generated eigenvalue distribution generated using L 1 ^ distribution requires rescaling 
through A — > and p(X) — » p(X)/(f> with 



r(i + M )cos(Y 



r(i 



(45) 



Below we show a sequence of results for p = 1.95, 1.75, 1.50, 1.25 and 1.00 with this rescaling. The comparison 
is for high statistics 400 x 400 samples (red). We have checked that the convergence is good already for 100 x 100 
samples, with no significant difference between N — 100, 200,400. The numerical results are also not sensitive to the 
choice (3 ^ 0. The agreement between the results following from the integral equations and the numerically generated 
spectra is perfect. This is true even for p, = 1, where in principle the arguments used in the derivation may not be 
valid. 



Numerical observation 

In the mean-field approximation 0] one can assume that there are no correlations between large eigenvalues of the 
Wigner-Levy random matrix. In this case the eigenvalue density takes the form: 

MX) = L C J?' m (\). (46) 

It is natural to ask how good this mean-field approximation is. This can be done by comparing the mean-field 
eigenvalue distribution p(X) (|4*6"|) to the eigenvalue distribution p(A) calculated by the inverse Hilbert transform ([77]) 
of the resolvent g(z) (|25|) as we did in previous section. We made this comparison numerically. The result of this 
numerical experiment was that within the numerical accuracy which we achieved the two curves representing p(X) and 
p(X) lied on top of each other in the whole studied range of A. Since our numerical codes are written in Mathematica 
we could push the numerical accuracy very far, being only limited by the execution time of the code. We have not 
seen any sign of deviation between the shapes of the two curves. This provides us with a strong numerical evidence 
that the mean-field argument Q gives an exact result but so far we have not managed to prove it. The value of the 
eigenvalue density p M (0) for A = can be calculated analytically for the mean-field density (|4*fJ|) . Rescaling the density 
Pfj,(X) — > p,j.{(j)\)/4> by the factor </> (|4"5"j) we eventually obtain 
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FIG. 6: The line represents the function Pn(0) (|47[) while the circles the values of p/i(0) computed numerically for p — 
1.0, 1.25, 1.5, 1.75, 1.95, 2.0 from the equation pTJl. 



We draw this function in Fig[6l In the same figure we also show points representing numerically evaluated values of 
the corresponding density /^(O) (|27p at some values of /i. Within the numerical accuracy p^(0) and p M (0) assume the 
same values. 

The physical meaning of the mean-field argument is that indeed one can think of the large eigenvalues as independent 
of each other. A similar observation has been made recently [l8j |. The mathematical meaning of the mean-field 
argument is more complex as we shall discuss below. 

Let us for brevity denote the function L^j^ \ x )> which is a function of two real arguments, by f(z,x) = 
L^S^ {x). The equation ((25)) can be now written as: 



/.:) = /- dx f -^^ (48) 



and (J27J) as 

p(X) = - 2 f dz-^-. (49) 

7T 2 Z-X 

When we insert (|48p into the last equation we have: 

^-^f dz f d \ z f{ X)(txy (50) 

A question is when this exact expression for p(A) is equal to the mean-field solution: p(X) = /(A, A). Recall the 
Poincare-Bertrand theorem. It tells us that the following equation holds 

/(A ' A) = * loo dz loo dx ( Z -x ){z -x) - jL dz ( Z -x )iz -xy (51) 



12 



We see that the density p(X) is given by the mean-field result: p(A) = p(A) = /(A, A) if the second term on the 
right-hand side of (p)Tj) vanishes. Unfortunately we have not managed to show that this is really the case for f(z,x) — 

(x) . One can however trivially observe that it would be the case if f(z,x) had the following form f(z,x) — 

f(x,x) = L^^'^ x '{x), and probably also if f(z,x) were a slowly varying function of z for z close to x, in which case 

the integral (|48[) would pick up only the contribution from f(x,x) leading to (|46[) . 

III. FREE RANDOM LEVY MATRICES 

Rotationally invariant measure 

Clearly Wigner Levy matrices are not rotationally invariant. In this section we shall discuss orthogonally (or 
unitary) invariant ensembles of Levy matrices. It can be shown that maximally random measures for such matrices 
have the form 0, [l6| : 

d» FR (H) = l[dH ij e- N ™W (52) 

i<j 

We shall be interested only in potentials with have tails which lead to eigenvalue distributions (spectral densities) 
with heavy tails p(X) ~ A~ 1_Al belonging to the Levy domain of attraction. A generic form of V(X) at asymptotic 
eigenvalues A is in this case 

V(X) =lnA 2 +C(l/A^) (53) 

In general the potential does not have to be an analytic function. We shall be interested here only in stable ensembles 
in the sense that the spectral measure (|52[) for the convolution of two independent and identical ensembles has the 
same form as the measure of the individual ensembles. In other words, the spectral measure for a matrix constructed 
as a sum of two independent matrices taken from the ensemble has exactly the same spectral measure (eigenvalue 
density) modulo linear transformations. 

It turns out that one can classify all the stable spectral measures thanks to the relation of the problem to free 
probability calculus. The matrix ensemble (|52p is in the large N limit a realization of free random variables 
so one can use theorems developed in free random probability [8]. In particular we can use the fact that in free 
probability theory stable laws are classified. They actually parallel stable laws of classical probability theory. In 
free probability the analogue of the logarithm of the characteristic function ([3]) is the R-transform, introduced by 
Voiculescu [12j. The R-transform linearizes the matrix convolution, generating spectral cumulants, which are additive 
under convolution. 

Stable laws in free probability 

The remarkable achievement by Bercovici and Voiculescu [8i] is an explicit derivation of all R-transforms defined 
by the equation R(G(z)) = z — 1/G(z) where G(z) is the resolvent for all free stable distributions. We just note that 
R(G(z)) is a sort of self-energy for rotationally symmetric FRL ensembles which is the analogue of (fl3|) which we 
previously defined for Wigner-Levy ensembles. It is self-averaging and additive. 

For stable laws R(z) is known. It can has either the trivial form R(z) = a or 

R(z) = bz^ 1 (54) 

where 0</i<2,6isa parameter which can be related to the stability index /i^the asymmetry parameter (3, and the 
range C known from the corresponding stable laws © of classical probability [1, [l| 

C e i{$-i)(M0)« for i< /i <2 
C e^+'iWM for < fj, < 1 ' 

In the marginal case: /it = 1, R(z) reads: 

R(z) = -iC(l+j3)-^—lnCz (55) 

The branch cut structure of R(z) is chosen in such a way that the upper complex half plane is mapped to itself. 
Recalling that R = z — 1/G in the large N limit, one finds that for the trivial case R(z) = 0, the resolvent: G(z) = z^ 1 
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FIG. 7: WL (black) versus FRL (red) for (j, = 1.95 



and the spectral distribution is a Dirac delta, p(X) = (5(A). Otherwise, on the upper half-plane, the resolvent fulfills 
an algebraic equation 



bG^{z) - zG{z) + 1 = 0, 
or in the marginal case {fi = 1): 

z+iC(l+0y\G{z) + —G(z) In CG(z) - 1 = 0. 



(56) 



(57) 

On the lower half-plane G(z) = G{z) [|. 

The equation for the resolvent (|56[) has explicit solutions only for the following values: p = 
1/4,1/3,1/2,2/3,3/4,4/3,3/2 and 2. In all other cases the equation is transcendental and one has to apply nu- 
merical procedures to unravel the spectral distribution. Again the form of the potential generating stable free Levy 
ensembles is highly non-trivial and is only known in few cases Q. We refer to Q for further references and discussions. 



Comparison of free Levy and Wigner-Levy spectra 

We present in Figures 7-11 several comparisons between the free random Levy spectra (FRL) following from the 
solution to the transcendental equation (red) and the random Levy spectra (BC) obtained by solving the coupled 
integral equations (black), for zero asymmetry ((3 = 0) and different tail indices fi. The FRL spectra are normalized 
to agree with the BC spectra in the tails of the distributions. We recall that the FRL spectra asymptote p(X) w 
sin(7r/i/2)/7r. The comparison in bulk shows that the spectra are similar, in particular close to the Gaussian limit 
/i = 2, where both approaches become equivalent. For smaller p there are differences. 

WL and FRL matrices represent two types of random matrices spectrally stable under matrix addition. For the 
WL matrices it follows from the measure, since each matrix elements is generated from a stable Levy distribution and 
therefore the sum of M WL matrices, scaled by is equivalent to the original WL ensemble. The important 

point is that the WL measure is not symmetric while the FRL one is. 



IV. SPECTRAL STABILITY AND MAXIMAL ENTROPY PRINCIPLE 



The matrix ensembles discussed in this paper are stable with respect to matrix addition in the sense that the 
eigenvalue distribution for the matrix constructed as a sum of two independent matrices from the original ensemble 
H = Hi + H2 is identical as the the original one up to a trivial rescaling. Wigner-Levy matrices are obviously stable, 
since the probability distribution for individual matrix elements of the sum Hij — Huj +l?2,-y is stable. A sum of two 
Wigner-Levy matrices is again a Wigner-Levy matrix. The Wigner-Levy matrices are not rotatationally invariant. 
This means in particular that the eigenvalue distribution itself does not provide the whole information about the 
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FIG. 11: WL (black) versus FRL (red) for /m = 1.00 



underlying matrix ensemble. Indeed, if O is a fixed orthogonal matrix, and H is a Wigner-Levy matrix, then the 
matrix OHO T is not anymore a Wigner-Levy matrix but it has exactly the same eigenvalue distribution as H . In 
other words, an ensemble of Wigner matrices is not maximally random among ensembles with the same eigenvalue 
distribution. One expects that maximally random ensemble with the given spectral properties should be rotationally 
invariant. In this case one also expects that the stability holds not only for the sum H = H\ + H2 but also for the sum 
of relatively rotated matrices: H = Hi + OH 2 T , where O is an arbitrary orthogonal matrix. It can be shown [l6| 
that the ensemble of random matrices which maximizes randomness (Shannon's entropy) for a given spectral density 
has the probability measure exactly of the form (|52|) as discussed here. 

Stable laws are important because they define domains of attractions. For example, if one thinks of a matrix 
addition one expects that a sum of many independent identically distributed random matrices H = Hi + • • • + H n 
should for n — > 00 become a random matrix from a stable ensemble. 

Maximally random spectrally stable ensembles which we discussed in the section on free random matrices play a 
special role since they can serve as an attraction point for the sums of iid rotationally invariant matrices. Moreover 
one expects that even for not rotationally invariant random matrices H, the sums of the form B = OiHiO\ + • ■ ■ + 
O n H n On where Oi are random orthogonal matrices, will for large n generate a maximally random matrix B from 
a spectrally stable ensemble. In this spirit one can expect that if ones adds many randomly rotated Wigner-Levy 
matrices: 



that for J\f — > 00 the matrices B should become rotationally invariant, maximally random with a distribution governed 
by the FRL symmetric distribution. In Figures 12-14 we show that this is indeed the case. The plots illustrate the 
two types of stability discussed above. In each case we generate N = 100 WL matrices and combine N — 100 of 
them either as a simple sum (black) or a rotated sum (red) with the appropriate scale factor. The plots represent the 
numerically measured spectra for the two cases. We present results for fi — 1.5, 1.25 and 1, which all show that a 
simple sum reproduces the BC result, while the rotated sum reproduces the symmetric FRL distribution. 



We have given a detailed analysis of the macroscopic limit of two distinct random matrix theories based on Levy 
type ensembles. The first one was put forward by Bouchaud and Cizeau [7j and uses a non-symmetric measure under 
the orthogonal group, and the second one was suggested by us [9j and uses a symmetric measure. 

After correcting the original analysis in Q, in particular our formulae (|23|) . (|2"Tj) replace (10b) and (12b) in Q, and 
their eq. (15) is replaced by the pair of "dispersion relations" (|22p and (|27p. we found perfect agreement between 
the analytical and numerical spectra. The WL measure is easy to implement numerically for arbitrary asymmetry 




(58) 



V. 



CONCLUSIONS 



A 

FIG. 12: WL (black) versus FRL (red) stability for jj, = 1.5 




A 

FIG. 14: WL (black) versus FRL (red) stability for fx = 1 
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parameter /? in the Levy distributions. The spectrum of WL matrices does not depend on j3 and remains symmetric 
and universal, depending only on /x. 

We have also shown that the spectra generated analytically for symmetric FRL matrices are similar to the ones 
generated from WL matrices. Unlike the WL ensemble, the FRL ensemble allows for both symmetric and asymmetric 
Levy distributions. Both ensembles are equally useful for addressing issues of recent interest [Til. Il5j]. 

Let us finish the paper with two remarks. 

• i. The application of free probability calculus to asymptotically free matrix realizations allows one to derive 
spectral density of the matrices from the underlying matrix ensemble but it does not tell one how to calculate 
eigenvalue correlation functions, or joint probabilities for many eigenvalues. Actually different matrix realiza- 
tions of free random variables may have completely different structure of eigenvalue correlations even if they are 
realizations of the same free random variables. To fix correlations or joint probabilities for two or more eigenval- 
ues, one has to introduce the concept of higher order freeness [Tjj. We think however that if one imposes on a 
matrix realization of free random variables an additional requirement that it has to be maximally random in the 
sense of maximizing Shannon's entropy [l6j then this aditional requirement automatically fixes the probability 
measure (j5"2")l for the ensemble and thus also all multi-eigenvalues correlations. 

• ii. Large eigenvalues behave differently for Wiener-Levy and maximally random free random Levy matrices 
discussed in this paper. As pointed out recently [18j ]. the largest eigenvalues fluctuate independently for Wigner- 
Levy ensemble, a little bit like in the mean-field argument [7[mentioned before, while for the maximally random 
matrix ensemble (|52|) even large eigenvalues are correlated [!| . 
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